using System;
using System.Text;
using System.Drawing;
using System.Drawing.Imaging;
using System.Windows.Forms ;
using System.Diagnostics ;
using System.Runtime.InteropServices;
using System.Collections;
using CommonNamespace;

namespace ROIImageProcessingNamespace
{
	/// <summary>
	/// Image processing routines using unmanaged code mostly
	/// because using setPixel and getPixel is SOOOO SLOOOOW
	/// </summary>
	/// 
	public unsafe class ImageProcessing
	{		
		//External Windows API
		[DllImport("kernel32.dll", EntryPoint="CopyMemory", CallingConvention=CallingConvention.StdCall, SetLastError=true)]
		public static extern void CopyMemory(void* PtrOut,void * PtrIn, int iByteCount);

		public enum ConvolutionKernelType
		{
			Blur = 0,
			GradientSobelHorizontal = 1,
			GradientSobelVertical = 2,
			GradientSobelHorizontalMagnitude = 3,
			GradientSobelVerticalMagnitude = 4,
			HighPass = 5,
			HighPassMagnitude = 6,
			Laplacian = 7,
			LaplacianMagnitude = 8
		}

		public struct ConvolutionKernel
		{
			public int [,] Matrix;
			public int Offset;
			public float NormalisationFactor;
			public bool TakeAbsoluteValue;
		}
		
		public enum BitmapMaskCombineOperation
		{
			Union = 0, 
			Intersect = 1,
			Xor = 2,
			Exclude = 3,
			Complement = 4
		}

		public enum BitmapMaskMorphologyOperation
		{
			Dilate = 0, 
			Erode = 1,
			Open = 2,
			Close = 3,
		}


		
		public static Bitmap Resize(Bitmap oBitmapIn, float fScale)
		{	
			int iHeight, iWidth;
			iHeight = (int) Math.Round(fScale * oBitmapIn.Height);
			iWidth = (int) Math.Round(fScale * oBitmapIn.Width );

			Bitmap oBitmapOut = new Bitmap(iWidth, iHeight);	
			Graphics g = Graphics.FromImage(oBitmapOut);	
			g.InterpolationMode = System.Drawing.Drawing2D.InterpolationMode.HighQualityBilinear;
			g.DrawImage(oBitmapIn, 0, 0, iWidth, iHeight);	
			return oBitmapOut;
		}
	
		public static Rectangle GetBitmapCenterRectangle(Bitmap objBitmap, int iWindowSize)
		{
			//Return a rectangle located in the center of the image. 
			Rectangle objRect = new System.Drawing.Rectangle((objBitmap.Width - iWindowSize) / 2,
				(objBitmap.Height - iWindowSize) / 2, 
				iWindowSize, iWindowSize);
			return(objRect);
		}
		
		
		public static unsafe int Average(Bitmap objBitmap, float[] sRGB)
		{
			//Get average pixel value in R*G*B*.
			//Return the number of pixels used.
			Rectangle objRect = new System.Drawing.Rectangle(0, 0, objBitmap.Width, objBitmap.Height);
			return Average(objBitmap, objRect, sRGB);
		}


		public static unsafe  int Average(Bitmap objBitmap, Rectangle objRect, float[] sRGB)
		{
			return Average(objBitmap, objRect, sRGB, null);
		}

		public static unsafe int Average(Bitmap objBitmap, float[] sRGB, Region objRegion)
		{
			Graphics g = Graphics.FromImage(objBitmap);
			
			if (! objRegion.IsEmpty(g))
			{
				RectangleF objRectF = objRegion.GetBounds(g);
				Rectangle objRect = new Rectangle((int) objRectF.X, (int)objRectF.Y, 
					(int)objRectF.Width, (int)objRectF.Height);
				return Average(objBitmap, objRect, sRGB, objRegion);	
			}
			else
			{
				return 0;
			}
		}	
			
		public static unsafe int Average(Bitmap objBitmap, Rectangle objRect, float[] sRGB, Region objRegion)
		{
			//Get average pixel value in R*G*B*.In Memory RGB values are stored in the order B-G-R!!!!!!
			//Getting a smaller rectangle with LockBits gets strange results, so we lock the whole image ...
			//Using a single to accumulate the averages can cause rounding errors, so we use an integer. 
			//This limits the use to images with about 8.5 megapixels, if they are all at 255...
			//Simply use long in that case ...
			int iPtrPixel0, iNrPixels = objRect.Height * objRect.Width, iRed = 0, iGreen = 0, iBlue = 0;
			Point Pt = new Point(0,0);

			if (iNrPixels > 0)
			{
				System.Drawing.Imaging.BitmapData bitmapData = objBitmap.LockBits(
					new Rectangle(0,0,objBitmap.Width, objBitmap.Height) , 
					System.Drawing.Imaging.ImageLockMode.ReadOnly ,
					System.Drawing.Imaging.PixelFormat.Format24bppRgb);																					
				iPtrPixel0 = bitmapData.Scan0.ToInt32();

				//Clip if necessary!
				objRect.Intersect(new Rectangle(0,0, objBitmap.Width, objBitmap.Height));
				
				//Byte pointer
				byte * byPixel = (byte *)(iPtrPixel0 + objRect.Top *  bitmapData.Stride + objRect.Left * 3);
				int iPtrJump = bitmapData.Stride - objRect.Width * 3;
				bool bUseRegion = false;
				if (objRegion != null) bUseRegion = true;
				
				//Recount the number of pixels, as the region may exclude a lot of them 
				iNrPixels = 0;
				for (Pt.Y = objRect.Top ; Pt.Y  < objRect.Bottom; Pt.Y  ++)
				{
					for (Pt.X  = objRect.Left; Pt.X  < objRect.Right;Pt.X ++)
					{
						if ((bUseRegion == false) || ((bUseRegion == true) && (objRegion.IsVisible(Pt))))
						{
							iBlue += *byPixel;
							byPixel ++;
							iGreen += *byPixel;
							byPixel ++;
							iRed += *byPixel;
							byPixel ++;	
							iNrPixels ++;
						}
						else
						{
							byPixel += 3;
						}
					}
					byPixel += iPtrJump;
				}

				objBitmap.UnlockBits(bitmapData);
				
				//Do computation of average in double, then convert to single
				sRGB[0] = (float)((double)iRed / iNrPixels);
				sRGB[1] = (float)((double)iGreen / iNrPixels);
				sRGB[2] = (float)((double)iBlue / iNrPixels);
			}
				
			return iNrPixels;				
		}

		public static unsafe void PolynomialTransform(Bitmap objBitmap, float[,] fLUT, float[,] fTransform)
		{
			//Do LUT and polynomial transform on whole image ...
			//Use unmanaged pointers to go FAST. In Memory RGB values are stored in the order B-G-R!!!!!!
			int x, y, iPtrPixel0;

			Rectangle rect = new Rectangle(0,0,objBitmap.Width,objBitmap.Height);
			System.Drawing.Imaging.BitmapData bitmapData = objBitmap.LockBits(rect, 
				ImageLockMode.ReadWrite,
				PixelFormat.Format24bppRgb);
																
			iPtrPixel0 = bitmapData.Scan0.ToInt32();

			//Byte pointer and corresponding managed array type
			byte * byPixel;
			int[] iRGB = new int[3];
			float[] fRGB = new float[3], fTransformRGB = new float[3];

			for (y = 0; y < bitmapData.Height; y ++)
			{
				byPixel = (byte *)(iPtrPixel0 + y * bitmapData.Stride);
				for (x = 0; x < bitmapData.Width; x ++)
				{
					iRGB[2] = *byPixel;
					iRGB[1] = *(byPixel + 1);
					iRGB[0] = *(byPixel + 2);

					//Apply LUT (R*G*B* to RGB)
					ColorUtil.RGBToLUTRGB(iRGB, fRGB, fLUT);

					//Apply transform
					ColorUtil.RGBTosRGB(fRGB, fTransformRGB, fTransform);

					//Apply Gamma correction
					ColorUtil.RGBToGammaRGB(fTransformRGB, iRGB);

					*byPixel = (byte) iRGB[2];
					byPixel ++;
					*byPixel = (byte) iRGB[1];
					byPixel ++;
					*byPixel = (byte) iRGB[0];
					byPixel ++;
				}
			}
			objBitmap.UnlockBits(bitmapData);
		}

		public static unsafe void Add(Bitmap objBitmap1, Bitmap objBitmap2)
		{
			Add(objBitmap1, objBitmap2, new Rectangle(0,0,objBitmap1.Width,objBitmap1.Height));
		}

		public static unsafe void Add(Bitmap objBitmap1, Bitmap objBitmap2, Rectangle rect)
		{
			//Add two images together, in R*G*B*. The result goes in the first image. Scale by 0.5.
			//Use unmanaged pointers to go FAST. In Memory RGB values are stored in the order B-G-R!!!!!!
			int x, y, iPtrPixel10, iPtrPixel20;

			System.Drawing.Imaging.BitmapData bitmapData1 = objBitmap1.LockBits(rect, 
				ImageLockMode.ReadWrite, PixelFormat.Format24bppRgb);
			System.Drawing.Imaging.BitmapData bitmapData2 = objBitmap2.LockBits(rect, 
				ImageLockMode.ReadOnly, PixelFormat.Format24bppRgb);

			iPtrPixel10 = bitmapData1.Scan0.ToInt32();
			iPtrPixel20 = bitmapData2.Scan0.ToInt32();

			//Byte pointer and corresponding managed array type
			byte * byPixel1, byPixel2;
			for (y = 0; y < bitmapData1.Height; y ++)
			{
				byPixel1 = (byte *)(iPtrPixel10 + y * bitmapData1.Stride);
				byPixel2 = (byte *)(iPtrPixel20 + y * bitmapData2.Stride);
				for (x = 0; x < bitmapData1.Width; x ++)
				{
					*byPixel1 = (byte) (*byPixel1 / 2 + *byPixel2 / 2);
					byPixel1 ++;
					byPixel2 ++;
					*byPixel1 = (byte) (*byPixel1 / 2 + *byPixel2 / 2);
					byPixel1 ++;
					byPixel2 ++;
					*byPixel1 = (byte) (*byPixel1 / 2 + *byPixel2 / 2);
					byPixel1 ++;
					byPixel2 ++;
				}
			}
			objBitmap1.UnlockBits(bitmapData1);
			objBitmap2.UnlockBits(bitmapData2);
		}

		public static ConvolutionKernel GetConvolutionKernel(ConvolutionKernelType KernelType, 
			int iKernelSize)
		{
			//A number of predefined kernels
			ConvolutionKernel objKernel = new ConvolutionKernel();

			//Allocate the kernel Matrix and Offset
			int i,j;

			switch (KernelType) 
			{
				case ConvolutionKernelType.Blur:
					objKernel.Matrix = new int[iKernelSize, iKernelSize];
					objKernel.NormalisationFactor = iKernelSize * iKernelSize;
					objKernel.TakeAbsoluteValue = false;
					for (i = 0; i < iKernelSize; i ++)
					{
						for (j = 0; j < iKernelSize; j ++)
						{
							objKernel.Matrix[i,j] = 1;
						}
					}
					break;
				case ConvolutionKernelType.GradientSobelHorizontal:
					objKernel.Matrix = new int[iKernelSize, iKernelSize];
					objKernel.NormalisationFactor = 6;
					objKernel.TakeAbsoluteValue = false;
					objKernel.Offset = 128;
					for (i = 0; i < iKernelSize; i ++)
					{
						objKernel.Matrix[i,0] = 1;
					}
					for (i = 0; i < iKernelSize; i ++)
					{
						objKernel.Matrix[i,iKernelSize-1] = -1;
					}
					break;
				case ConvolutionKernelType.GradientSobelVertical:
					objKernel.Matrix = new int[iKernelSize, iKernelSize];
					objKernel.NormalisationFactor = 6;
					objKernel.TakeAbsoluteValue = false;
					objKernel.Offset = 128;
					for (j = 0; j < iKernelSize; j ++)
					{
						objKernel.Matrix[0,j] = 1;
					}
					for (j = 0; j < iKernelSize; j ++)
					{
						objKernel.Matrix[iKernelSize-1, j] = -1;
					}
					break;
				case ConvolutionKernelType.GradientSobelHorizontalMagnitude:
					objKernel.Matrix = new int[iKernelSize, iKernelSize];
					objKernel.NormalisationFactor = 3;
					objKernel.TakeAbsoluteValue = true;
					for (i = 0; i < iKernelSize; i ++)
					{
						objKernel.Matrix[i,0] = 1;
					}
					for (i = 0; i < iKernelSize; i ++)
					{
						objKernel.Matrix[i,iKernelSize-1] = -1;
					}
					break;
				case ConvolutionKernelType.GradientSobelVerticalMagnitude:
					objKernel.Matrix = new int[iKernelSize, iKernelSize];
					objKernel.NormalisationFactor = 3;
					objKernel.TakeAbsoluteValue = true;
					for (j = 0; j < iKernelSize; j ++)
					{
						objKernel.Matrix[0,j] = 1;
					}
					for (j = 0; j < iKernelSize; j ++)
					{
						objKernel.Matrix[iKernelSize-1, j] = -1;
					}
					break;
				case ConvolutionKernelType.HighPass:
					objKernel.TakeAbsoluteValue = false;
					objKernel.Offset = 128;
				switch(iKernelSize)
				{
					case 3:
						objKernel.Matrix = new int[,] {{-1,-1,-1},{-1,9,-1},{-1,-1,-1}};
						objKernel.NormalisationFactor = 18;
						break;
					case 5:
						objKernel.Matrix = new int[,] {{0,-1,-1,-1,0},{-1,2,-4,2,-1},{1,-4,13,-4,1}, {0,-1,-1,-1,0},{-1,2,-4,2,-1}};
						objKernel.NormalisationFactor = 46;
						break;
					default:
						break;
				}
					break;
				case ConvolutionKernelType.HighPassMagnitude :
					objKernel.Offset = 0;
					objKernel.TakeAbsoluteValue = true;
				switch(iKernelSize)
				{
					case 3:
						objKernel.Matrix = new int[,] {{-1,-1,-1},{-1,9,-1},{-1,-1,-1}};
						objKernel.NormalisationFactor = 9;
						break;
					case 5:
						objKernel.Matrix = new int[,] {{0,-1,-1,-1,0},{-1,2,-4,2,-1},{1,-4,13,-4,1}, {0,-1,-1,-1,0},{-1,2,-4,2,-1}};
						objKernel.NormalisationFactor = 23;
						break;
					default:
						break;
				}
					break;
				case ConvolutionKernelType.Laplacian :
					objKernel.Offset = 128;
					objKernel.TakeAbsoluteValue = false;
				switch(iKernelSize)
				{
					case 3:
						objKernel.Matrix = new int[,] {{0,-1,0},{-1,4,-1},{0,-1,0}};
						objKernel.NormalisationFactor = 8;
						break;
					case 5:
						objKernel.Matrix = new int[,] {{-1,-1,-1,-1,-1},{-1,-1,-1,-1,-1},{-1,-1,24,-1,-1}, {-1,-1,-1,-1,-1},{-1,-1,-1,-1,-1}};
						objKernel.NormalisationFactor = 48;
						break;
					default:
						break;
				}
					break;
				case ConvolutionKernelType.LaplacianMagnitude :
					objKernel.Offset = 0;
					objKernel.TakeAbsoluteValue = true;
				switch(iKernelSize)
				{
					case 3:
						objKernel.Matrix = new int[,] {{0,-1,0},{-1,4,-1},{0,-1,0}};
						objKernel.NormalisationFactor = 4;
						break;
					case 5:
						objKernel.Matrix = new int[,] {{-1,-1,-1,-1,-1},{-1,-1,-1,-1,-1},{-1,-1,24,-1,-1}, {-1,-1,-1,-1,-1},{-1,-1,-1,-1,-1}};
						objKernel.NormalisationFactor = 24;
						break;
					default:
						break;
				}
					break;
				default:
					break;
			}
			return objKernel;
		}

		public static unsafe void Convolve(Bitmap objBitmap, 
			ConvolutionKernelType KernelType, 
			int iKernelSize)
		{
			//No output bitmap specified, make a copy of input bitmap, and return result in that
			//If repeated operations on same image, it might be faster to allocate output bitmap
			//in calling application
			Bitmap objSrcBitmap = new Bitmap(objBitmap);
			Convolve(objSrcBitmap, GetConvolutionKernel(KernelType, iKernelSize), objBitmap, 
				new Rectangle(0,0,objBitmap.Width,objBitmap.Height));
		}

		public static unsafe void Convolve(Bitmap objBitmap, 
			ConvolutionKernel objKernel)
		{
			//No output bitmap specified, make a copy of input bitmap, and return result in that
			//If repeated operations on same image, it might be faster to allocate output bitmap
			//in calling application
			Bitmap objSrcBitmap = new Bitmap(objBitmap);
			Convolve(objSrcBitmap, objKernel, objBitmap, new Rectangle(0,0,objBitmap.Width,objBitmap.Height));
		}

		public static unsafe void Convolve(Bitmap objBitmap, 
			ConvolutionKernelType KernelType, 
			int iKernelSize, 
			Bitmap objDstBitmap)
		{
			Convolve(objBitmap, GetConvolutionKernel(KernelType, iKernelSize), objDstBitmap, 
				new Rectangle(0,0,objBitmap.Width,objBitmap.Height));
		}

		public static unsafe void Convolve(Bitmap objSrcBitmap, 
			ConvolutionKernel objKernel, 
			Bitmap objDstBitmap)
		{
			Convolve(objSrcBitmap, objKernel, objDstBitmap, new Rectangle(0,0,objSrcBitmap.Width,objSrcBitmap.Height));
		}

		public static unsafe void Convolve(Bitmap objBitmap, 
			ConvolutionKernelType KernelType, 
			int iKernelSize,
			Rectangle rect)
		{
			//No output bitmap specified, make a copy of input bitmap, and return result in that
			//If repeated operations on same image, it might be faster to allocate output bitmap
			//in calling application
			Bitmap objSrcBitmap = new Bitmap(objBitmap);
			Convolve(objSrcBitmap, GetConvolutionKernel(KernelType, iKernelSize), objBitmap, rect);
		}

		public static unsafe void Convolve(Bitmap objBitmap, 
			ConvolutionKernel objKernel, 
			Rectangle rect)
		{
			//No output bitmap specified, make a copy of input bitmap, and return result in that
			//If repeated operations on same image, it might be faster to allocate output bitmap
			//in calling application
			Bitmap objSrcBitmap = new Bitmap(objBitmap);
			Convolve(objSrcBitmap, objKernel, objBitmap, rect);
		}

		public static unsafe void Convolve(Bitmap objSrcBitmap, 
			ConvolutionKernelType KernelType, 
			int iKernelSize, 
			Bitmap objDstBitmap,
			Rectangle rect)
		{
			Convolve(objSrcBitmap, GetConvolutionKernel(KernelType, iKernelSize), objDstBitmap, rect);
		}


		public static unsafe void Convolve(Bitmap objSrcBitmap, 
			ConvolutionKernel objKernel, 
			Bitmap objDstBitmap, 
			Rectangle rect)
		{
			//Do convolution. Output bitmap must be allocated and the same size as input bitmap
			//Output bitmap MUST be different as input bitmap ...
			//Convolution kernel must be odd-sized (3x3, 5x5, 7x7 etc...), and NormalisationFactorALIZED!
			//Use unmanaged pointers to go FAST. In Memory RGB values are stored in the order B-G-R!!!!!!

			System.Drawing.Imaging.BitmapData bitmapSrcData = objSrcBitmap.LockBits(rect, 
				System.Drawing.Imaging.ImageLockMode.ReadOnly ,
				System.Drawing.Imaging.PixelFormat.Format24bppRgb);
			System.Drawing.Imaging.BitmapData bitmapDstData = objDstBitmap.LockBits(rect, 
				System.Drawing.Imaging.ImageLockMode.WriteOnly,
				System.Drawing.Imaging.PixelFormat.Format24bppRgb);
																					
			int iPtrSrcPixel0 = bitmapSrcData.Scan0.ToInt32();
			int iPtrDstPixel0 = bitmapDstData.Scan0.ToInt32();

			//Byte pointer and corresponding managed array type
			byte * bySrcPixel, byDstPixel;
			float[] fSafePixel = new float[3];
			int x, y, i, j, iKernelHalfSize = objKernel.Matrix.GetUpperBound(0) / 2;
			int iKernelSize = 2 * iKernelHalfSize + 1;

			for (y = 0; y < bitmapSrcData.Height; y ++)
			{
				byDstPixel = (byte *)(iPtrDstPixel0 + y * bitmapDstData.Stride);
				for (x = 0; x < bitmapSrcData.Width; x ++)
				{
					fSafePixel[2] = 0;
					fSafePixel[1] = 0;
					fSafePixel[0] = 0;
					for (i = 0; i < iKernelSize; i ++)
					{
						for (j = 0; j < iKernelSize; j ++)
						{
							int x2, y2;
							x2 = x + j;
							y2 = y + i;
							if (x2 < iKernelHalfSize) x2 = iKernelHalfSize;
							if (x2 > (bitmapSrcData.Width - iKernelHalfSize)) x2 = bitmapSrcData.Width - iKernelHalfSize;
							if (y2 < iKernelHalfSize) y2 = iKernelHalfSize;
							if (y2 > (bitmapSrcData.Height - iKernelHalfSize)) y2 = bitmapSrcData.Height - iKernelHalfSize;
							
							bySrcPixel = (byte *)(iPtrSrcPixel0 + (y2 - iKernelHalfSize) * bitmapSrcData.Stride + 
								(x2 - iKernelHalfSize) * 3);

							fSafePixel[2] += *bySrcPixel * objKernel.Matrix[i,j];
							fSafePixel[1] += *(bySrcPixel + 1) * objKernel.Matrix[i,j];
							fSafePixel[0] += *(bySrcPixel + 2) * objKernel.Matrix[i,j];
						}
					}

					if (objKernel.TakeAbsoluteValue == true)
					{
						*byDstPixel = (byte) Math.Round(Math.Abs(fSafePixel[2] / objKernel.NormalisationFactor + objKernel.Offset));
						byDstPixel ++;
						*byDstPixel = (byte) Math.Round(Math.Abs(fSafePixel[1] / objKernel.NormalisationFactor + objKernel.Offset));
						byDstPixel ++;
						*byDstPixel = (byte) Math.Round(Math.Abs(fSafePixel[0] / objKernel.NormalisationFactor + objKernel.Offset));
						byDstPixel ++;
					}
					else
					{
						*byDstPixel = (byte) Math.Round(fSafePixel[2] / objKernel.NormalisationFactor + objKernel.Offset);
						byDstPixel ++;
						*byDstPixel = (byte) Math.Round(fSafePixel[1] / objKernel.NormalisationFactor + objKernel.Offset);
						byDstPixel ++;
						*byDstPixel = (byte) Math.Round(fSafePixel[0] / objKernel.NormalisationFactor + objKernel.Offset);
						byDstPixel ++;
					}
				}
			}
			objDstBitmap.UnlockBits(bitmapDstData);
			objSrcBitmap.UnlockBits(bitmapSrcData);
		}


		public static unsafe void MedianFilter(Bitmap objBitmap, int iKernelSize)
		{
			//No output bitmap specified, make a copy of input bitmap, and return result in that
			//If repeated operations on same image, it might be faster to allocate output bitmap
			//in calling application
			Bitmap objSrcBitmap = new Bitmap(objBitmap);
			MedianFilter(objSrcBitmap, iKernelSize, objBitmap, new Rectangle(0,0,objBitmap.Width,objBitmap.Height));
		}

		public static unsafe void MedianFilter(Bitmap objSrcBitmap, int iKernelSize, Bitmap objDstBitmap)
		{
			MedianFilter(objSrcBitmap, iKernelSize, objDstBitmap, new Rectangle(0,0,objSrcBitmap.Width,objSrcBitmap.Height));
		}

		public static unsafe void MedianFilter(Bitmap objBitmap, int iKernelSize, Rectangle rect)
		{
			//No output bitmap specified, make a copy of input bitmap, and return result in that
			//If repeated operations on same image, it might be faster to allocate output bitmap
			//in calling application
			Bitmap objSrcBitmap = new Bitmap(objBitmap);
			MedianFilter(objSrcBitmap, iKernelSize, objBitmap, rect);
		}

		public static unsafe void MedianFilter(Bitmap objSrcBitmap, int iKernelSize, Bitmap objDstBitmap, Rectangle rect)
		{
			//Do median filtering. Output bitmap must be allocated and the same size as input bitmap
			//Output bitmap MUST be different as input bitmap ...
			//Use unmanaged pointers to go FAST. In Memory RGB values are stored in the order B-G-R!!!!!!
			System.Drawing.Imaging.BitmapData bitmapSrcData = objSrcBitmap.LockBits(rect, 
				System.Drawing.Imaging.ImageLockMode.ReadOnly ,
				System.Drawing.Imaging.PixelFormat.Format24bppRgb);
			System.Drawing.Imaging.BitmapData bitmapDstData = objDstBitmap.LockBits(rect, 
				System.Drawing.Imaging.ImageLockMode.WriteOnly,
				System.Drawing.Imaging.PixelFormat.Format24bppRgb);
																					
			int iPtrSrcPixel0 = bitmapSrcData.Scan0.ToInt32();
			int iPtrDstPixel0 = bitmapDstData.Scan0.ToInt32();

			//Byte pointer and corresponding managed array type
			byte * bySrcPixel, byDstPixel;
			int x, y, i, j, iKernelHalfSize = (iKernelSize - 1) / 2, iKernelPixelNr;
			int iNrKernelPixels = iKernelSize * iKernelSize, iKernelCentralPixel = (iNrKernelPixels - 1) / 2;
			float[] fSafePixel = new float[2];
			float[] fSafePixelRed = new float[iNrKernelPixels];
			float[] fSafePixelGreen = new float[iNrKernelPixels];
			float[] fSafePixelBlue = new float[iNrKernelPixels];

			for (y = 0; y < bitmapSrcData.Height; y ++)
			{
				byDstPixel = (byte *)(iPtrDstPixel0 + y * bitmapDstData.Stride);
				for (x = 0; x < bitmapSrcData.Width; x ++)
				{
					iKernelPixelNr = 0;
					for (i = 0; i < iKernelSize; i ++)
					{
						for (j = 0; j < iKernelSize; j ++)
						{
							int x2, y2;
							x2 = x + j;
							y2 = y + i;
							if (x2 < iKernelHalfSize) x2 = iKernelHalfSize;
							if (x2 > (bitmapSrcData.Width - iKernelHalfSize)) x2 = bitmapSrcData.Width - iKernelHalfSize;
							if (y2 < iKernelHalfSize) y2 = iKernelHalfSize;
							if (y2 > (bitmapSrcData.Height - iKernelHalfSize)) y2 = bitmapSrcData.Height - iKernelHalfSize;
							
							bySrcPixel = (byte *)(iPtrSrcPixel0 + (y2 - iKernelHalfSize) * bitmapSrcData.Stride + 
								(x2 - iKernelHalfSize) * 3);

							fSafePixelBlue[iKernelPixelNr] = *bySrcPixel;
							fSafePixelGreen[iKernelPixelNr] = *(bySrcPixel + 1);
							fSafePixelRed[iKernelPixelNr] = *(bySrcPixel + 2);
							
							iKernelPixelNr ++;
						}
					}
					
					//Sort the arrays
					Array.Sort(fSafePixelBlue);
					Array.Sort(fSafePixelRed);
					Array.Sort(fSafePixelGreen);
					
					*byDstPixel = (byte) (fSafePixelBlue[iKernelCentralPixel]);
					byDstPixel ++;
					*byDstPixel = (byte) (fSafePixelGreen[iKernelCentralPixel]);
					byDstPixel ++;
					*byDstPixel = (byte) (fSafePixelRed[iKernelCentralPixel]);
					byDstPixel ++;
				}
			}
			objDstBitmap.UnlockBits(bitmapDstData);
			objSrcBitmap.UnlockBits(bitmapSrcData);
		}

		public static unsafe void Multiply(Bitmap objBitmap, float[] fFactor)
		{
			//Multiply image data with fFactor
			//Use unmanaged pointers to go FAST. In Memory RGB values are stored in the order B-G-R!!!!!!
			Rectangle rect = new Rectangle(0,0,objBitmap.Width, objBitmap.Height);
			System.Drawing.Imaging.BitmapData bitmapData = objBitmap.LockBits(rect, 
				System.Drawing.Imaging.ImageLockMode.ReadOnly ,
				System.Drawing.Imaging.PixelFormat.Format24bppRgb);
																					
			int iPtrPixel0 = bitmapData.Scan0.ToInt32();

			//Byte pointer and corresponding managed array type
			byte * byPixel;
			int x, y;
			for (y = 0; y < bitmapData.Height; y ++)
			{
				byPixel = (byte *)(iPtrPixel0 + y * bitmapData.Stride);
				for (x = 0; x < bitmapData.Width; x ++)
				{
					*byPixel = (byte) (*byPixel * fFactor[2]);
					byPixel ++;
					*byPixel = (byte) (*byPixel * fFactor[1]);
					byPixel ++;
					*byPixel = (byte) (*byPixel * fFactor[0]);
					byPixel ++;
				}
			}
			objBitmap.UnlockBits(bitmapData);
		}


		public static unsafe void Extrema(Bitmap objBitmap, byte[] byMinRGB, byte[] byMaxRGB)
		{
			Extrema(objBitmap, byMinRGB, byMaxRGB, new Rectangle(0,0,objBitmap.Width,objBitmap.Height));
		}

		public static unsafe void Extrema(Bitmap objBitmap, byte[] byMinRGB, byte[] byMaxRGB, Rectangle rect)
		{
			//Get minimal/maximal  R*G*B* value of image
			//Use unmanaged pointers to go FAST. In Memory RGB values are stored in the order B-G-R!!!!!!

			System.Drawing.Imaging.BitmapData bitmapData = objBitmap.LockBits(rect, 
				System.Drawing.Imaging.ImageLockMode.ReadOnly ,
				System.Drawing.Imaging.PixelFormat.Format24bppRgb);
																					
			int iPtrPixel0 = bitmapData.Scan0.ToInt32();

			//Byte pointer and corresponding managed array type
			byte * byPixel;
			int x, y;
			byMinRGB[0] = 255;
			byMinRGB[1] = 255;
			byMinRGB[2] = 255;
			byMaxRGB[0] = 0;
			byMaxRGB[1] = 0;
			byMaxRGB[2] = 0;
			for (y = 0; y < bitmapData.Height; y ++)
			{
				byPixel = (byte *)(iPtrPixel0 + y * bitmapData.Stride);
				for (x = 0; x < bitmapData.Width; x ++)
				{
					if (*byPixel < byMinRGB[2]) byMinRGB[2] = (byte) *byPixel;
					if (*byPixel > byMaxRGB[2]) byMaxRGB[2] = (byte) *byPixel;
					byPixel ++;
					if (*byPixel < byMinRGB[1]) byMinRGB[1] = (byte) *byPixel;
					if (*byPixel > byMaxRGB[1]) byMaxRGB[1] = (byte) *byPixel;
					byPixel ++;
					if (*byPixel < byMinRGB[0]) byMinRGB[0] = (byte) *byPixel;
					if (*byPixel > byMaxRGB[0]) byMaxRGB[0] = (byte) *byPixel;
					byPixel ++;
				}
			}
			objBitmap.UnlockBits(bitmapData);
		}

  #region Region functions
		//These region functions might well be obsolete, as it is MUCH faster to work with masks (1-bit images)
		public static bool IsConnected(float X, float Y, Region oRegion, Bitmap objBitmap)
		{
			//Return if a point from a region is 4-connected within that region
			//A point on the border of the image will only be tested for 3 or 2 connectedness.
			//This is a computationally expensive function!!!
			float X1,X2,Y1,Y2;

			X1 = Math.Max(0, X - 1);
			Y1 = Math.Max(0, Y - 1);
			X2 = Math.Min(objBitmap.Width - 1, X + 1);
			Y2 = Math.Min(objBitmap.Height - 1, Y + 1);

			if (oRegion.IsVisible(X1, Y) && oRegion.IsVisible(X2, Y) && 
				oRegion.IsVisible(X, Y1) && oRegion.IsVisible(X, Y2))
			{
				return true;
			}
			else
			{
				return false;
			}
		}

		private static void Erode(Bitmap objBitmap, Region objRegion)
		{
			Erode(objBitmap, objRegion, new RectangleF(0,0,objBitmap.Width, objBitmap.Height));
		}

		private static void Erode(Bitmap objBitmap, Region objRegion, RectangleF objRect)
		{
			//Erode (and dilate) a region, using 4-connection. We do not erode the border of the image!!!
			//This is a computationally very expensive function, so we try to optimize it as much as possible
			//Most expensive are the Region.IsVisible calls.
			//Debug.WriteLine("Eroding region in rectangle " + objRect.ToString(), "ErodeRegion");
			RectangleF objExcludeRect = new Rectangle(0,0,1,1);
			Region oOriginalRgn = objRegion.Clone();
 
			int x, y, xMin, yMin, xMax, yMax, iNrRemovedPixels = 0;
			yMin = (int) Math.Max(objRect.Top,1);
			yMax = (int) Math.Min(objRect.Bottom, objBitmap.Height - 2);
			xMin = (int) Math.Max(objRect.Left,1);
			xMax = (int) Math.Min(objRect.Right, objBitmap.Width - 2);

			//The algoritm will store the visibility of the next line to avoid retesting the same location 
			//more than once.
			bool[,] bPtVisible = new bool[3,xMax - xMin + 2];
			int iCurrentLineIndex = 1, iNextLineIndex = 2, iPreviousLineIndex = 0;

			y = yMin;
			for (x = xMin - 1; x <= xMax; x ++)
			{
				bPtVisible[iPreviousLineIndex, x - xMin + 1] = oOriginalRgn.IsVisible((float)x,(float)y - 1);
				bPtVisible[iCurrentLineIndex, x - xMin + 1] = oOriginalRgn.IsVisible((float)x,(float)y);
			}

			for (y = yMin; y < yMax; y ++)
			{
				//Store first and last point of of next line
				bPtVisible[iNextLineIndex, 0] = oOriginalRgn.IsVisible((float)0,(float)y + 1);
				bPtVisible[iNextLineIndex, xMax - xMin + 1] = oOriginalRgn.IsVisible((float)xMax,(float)y + 1);

				for (x = xMin; x < xMax; x ++)
				{
					int iCurrentPtIndex = x - xMin + 1;

					//Store point of next line
					bPtVisible[iNextLineIndex, iCurrentPtIndex] = oOriginalRgn.IsVisible((float)x,(float)y + 1);

					//Test current point if required
					if (bPtVisible[iCurrentLineIndex, iCurrentPtIndex])
					{
						//Point is in region, test for connectedness
						if (bPtVisible[iPreviousLineIndex, iCurrentPtIndex] && bPtVisible[iCurrentLineIndex, iCurrentPtIndex - 1] &&
							bPtVisible[iCurrentLineIndex, iCurrentPtIndex + 1] && bPtVisible[iNextLineIndex, iCurrentPtIndex])
						{
							//Connected, do not remove
						}
						else
						{
							//Remove point from region
							//Debug.WriteLine("Found border point in region at " + x + "," + y, "ErodeRegion");
							objExcludeRect.X = x;
							objExcludeRect.Y = y;
							objRegion.Exclude(objExcludeRect);
							iNrRemovedPixels ++;
						}
					}
				}

				//Juggle indices ...
				switch (iCurrentLineIndex)
				{         
					case 0:   
						iCurrentLineIndex = 1;
						iPreviousLineIndex = 0;
						iNextLineIndex = 2;
						break;
					case 1:
						iCurrentLineIndex = 2;
						iPreviousLineIndex = 1;
						iNextLineIndex = 0;
						break;
					case 2:
						iCurrentLineIndex = 0;
						iPreviousLineIndex = 2;
						iNextLineIndex = 1;
						break;
				}

			}

			//Dispose of region
			oOriginalRgn.Dispose();
			//Debug.WriteLine("Removed " + iNrRemovedPixels + " pixels", "ErodeRegion");
		}

		public static void Dilate(Bitmap objBitmap,Region oRegion)
		{
			//Dilate a region, using 4-connection. Can be implemented by erosion of the inverse region
			//When doing dilation using the erode routine, we can limit the 
			//bounding box enormously because we know the location of the region we inverted.
			Graphics g = Graphics.FromImage(objBitmap);
			RectangleF oRectOrignalRegion = oRegion.GetBounds(g);
			oRectOrignalRegion.Inflate(2,2);

			RectangleF oRect = new RectangleF(0,0,objBitmap.Width, objBitmap.Height);
			//Debug.WriteLine("Dilating region in rectangle " + oRect.ToString(), "DilateRegion (using ErodeRegion on complement)");
			oRegion.Complement(oRect);
			Erode(objBitmap, oRegion, oRectOrignalRegion);
			oRegion.Complement(oRect);
			g.Dispose();
		}

		public static void Close(Bitmap oBitmap,Region oRegion)
		{
			//Close is dilation followed by erosion, and close small gaps and smooth the contour
			Dilate(oBitmap, oRegion);
			Erode(oBitmap, oRegion);
		}

		public static void Open(Bitmap oBitmap,Region oRegion)
		{
			//Open is erosion followed by dilation, and removes small objects and smooth the contour
			Erode(oBitmap, oRegion);
			Dilate(oBitmap, oRegion);
		}

		public static Region GetBorderRegion(Region objRegion, Rectangle objRect) 
		{
			Region oBorderRegion = new Region();
			GetBorderRegion(objRegion, oBorderRegion, objRect);
			return (oBorderRegion);
		}

		public static float MeasureRegion(Bitmap objBitmap, Region oRegion, float [] sLab)
		{
			//Return color difference with given CIE Lab value
			float [] sResLab = new float[3];
			float [] sRGB = new float[3];
			MeasureRegion(objBitmap, oRegion, sRGB, sResLab);
			return (ColorUtil.CIELabTodE(sLab, sResLab));
		}

		public static void MeasureRegion(Bitmap objBitmap, Region oRegion, float [] sGammaRGB, float [] sLab)
		{
			float sA = 0;
			MeasureRegion(objBitmap, oRegion, sGammaRGB, sLab, ref sA);
		}

		public static void MeasureRegion(Bitmap objBitmap, Region oRegion, float [] sGammaRGB, float [] sLab, ref float sArea)
		{
			//Return area in cm^2 and average color of region, both sRGB and CIE Lab
			int iNrPixels = Average(objBitmap, sGammaRGB, oRegion);
			sArea = iNrPixels * (float) Math.Pow((2.54 / objBitmap.HorizontalResolution), 2);
			ColorUtil oColorUtil = new ColorUtil();
			float [] sRGB = new float[3];
			ColorUtil.GammaRGBToRGB(sGammaRGB, sRGB);
			oColorUtil.sRGBToCIELab(sRGB, sLab);
		}

		public static void GetBorderRegion(Region objRegion, Region oBorderRegion, Rectangle objRect)
		{
			//Determine the border region of a region
			//This is a computationally very expensive function, so we try to optimize it as much as possible
			//Debug.WriteLine("Eroding region in rectangle " + objRect.ToString(), "ErodeRegion");
			RectangleF objIncludeRect = new Rectangle(0,0,1,1);
			Region oOriginalRgn = objRegion.Clone();
 
			int x, y, xMin, yMin, xMax, yMax, iNrBorderPixels = 0;
			yMin = (int) objRect.Top + 1;
			yMax = (int) objRect.Bottom - 1;
			xMin = (int) objRect.Left + 1;
			xMax = (int) objRect.Right - 1;

			//The algoritm will store the visibility of the next line to avoid retesting the same location 
			//more than once.
			bool[,] bPtVisible = new bool[3,xMax - xMin + 2];
			int iCurrentLineIndex = 1, iNextLineIndex = 2, iPreviousLineIndex = 0;

			y = yMin;
			for (x = xMin - 1; x <= xMax; x ++)
			{
				bPtVisible[iPreviousLineIndex, x - xMin + 1] = oOriginalRgn.IsVisible((float)x,(float)y - 1);
				bPtVisible[iCurrentLineIndex, x - xMin + 1] = oOriginalRgn.IsVisible((float)x,(float)y);
			}

			for (y = yMin; y < yMax; y ++)
			{
				//Store first and last point of of next line
				bPtVisible[iNextLineIndex, 0] = oOriginalRgn.IsVisible((float)0,(float)y + 1);
				bPtVisible[iNextLineIndex, xMax - xMin + 1] = oOriginalRgn.IsVisible((float)xMax,(float)y + 1);

				for (x = xMin; x < xMax; x ++)
				{
					int iCurrentPtIndex = x - xMin + 1;

					//Store point of next line
					bPtVisible[iNextLineIndex, iCurrentPtIndex] = oOriginalRgn.IsVisible((float)x,(float)y + 1);

					//Test current point if required
					if (bPtVisible[iCurrentLineIndex, iCurrentPtIndex])
					{
						//Point is in region, test for connectedness
						if (!bPtVisible[iPreviousLineIndex, iCurrentPtIndex] || !bPtVisible[iCurrentLineIndex, iCurrentPtIndex - 1] ||
							!bPtVisible[iCurrentLineIndex, iCurrentPtIndex + 1] || !bPtVisible[iNextLineIndex, iCurrentPtIndex])
						{
							//Debug.WriteLine("Found border point in region at " + x + "," + y, "ErodeRegion");
							objIncludeRect.X = x;
							objIncludeRect.Y = y;
							objRegion.Union(objIncludeRect);
							iNrBorderPixels ++;
						}
					}
				}

				//Juggle indices ...
				switch (iCurrentLineIndex)
				{         
					case 0:   
						iCurrentLineIndex = 1;
						iPreviousLineIndex = 0;
						iNextLineIndex = 2;
						break;
					case 1:
						iCurrentLineIndex = 2;
						iPreviousLineIndex = 1;
						iNextLineIndex = 0;
						break;
					case 2:
						iCurrentLineIndex = 0;
						iPreviousLineIndex = 2;
						iNextLineIndex = 1;
						break;
				}

			}
			Debug.WriteLine("Found " + iNrBorderPixels + " border pixels", "GetBorderRegion");
		}

		public static bool Grow(Bitmap objBitmap, Region oRegion, float [] sSeedLab, float sdEThreshold, int iMaxNrIterations)
		{
			//Grow a region using given ColorUtil and dE threshold given.
			//There is not much point in using unmanaged code as pixel access will be random
			//Returns after given number of iterations, with true when finished, false otherwise
			//This allows update of UI.

			//Repeated dilation is computationally very expensive, so 
			//we try several tricks to speed things up ...

			Region oBorderRegion;
			RectangleF objPtRect = new RectangleF(0,0,1,1), objRect = new RectangleF(0,0,1,1);
			Color oColor;
			ColorUtil oColorClass = new ColorUtil();
			float[] sRGB = new float[3];
			float[] sLab = new float[3];
			float fdE;
			Graphics g = Graphics.FromImage(objBitmap);
			int iIteration = 0, iNewPoints = 0, iIterationNewPoints;
	
			do
			{
				//Determine exterior border of region.
				oBorderRegion = oRegion.Clone();
				Dilate(objBitmap, oBorderRegion);
				oBorderRegion.Exclude(oRegion);
				
				//Go through all the pixels, and add them to the original region if their ColorUtil 
				//is close enough to that of the seed point
				if (iIteration % 16 == 0)
				{
					objRect = oBorderRegion.GetBounds(g);
				}
				else
				{
					objRect.Inflate(1,1);
				}

				Debug.WriteLine("Checking border pixels in rectangle " + objRect.ToString() + " for iteration " + iIteration , "GrowRegion");
 
				int x, y, xMin, yMin, xMax, yMax, iNrBorderPixels = 0, iNrCheckedPixels = 0;
				yMin = (int) objRect.Top;
				yMax = (int) objRect.Bottom;
				xMin = (int) objRect.Left;
				xMax = (int) objRect.Right;
 
				iIterationNewPoints = 0;
				for (y = yMin; y < yMax; y ++)
				{
					for (x = xMin; x < xMax; x ++)
					{
						if (oBorderRegion.IsVisible((float)x,(float)y))
						{
							//Check ColorUtil
							oColor = objBitmap.GetPixel(x,y);
							ColorUtil.ColorToRGB(oColor, sRGB);
							oColorClass.sRGBToCIELab(sRGB, sLab);
		
							fdE = ColorUtil.CIELabTodE(sLab, sSeedLab);
						
							//Debug.WriteLine("Pixel at " + x + "," + y + " has dE " + fdE + " with seed", "GrowRegion");

							//Add to region if below threshold
							if (fdE < sdEThreshold)
							{
								objPtRect.X = (float) (x - (float)0.5);
								objPtRect.Y = (float) (y - (float)0.5);
								oRegion.Union(objPtRect);
								iIterationNewPoints ++;
							}

							iNrBorderPixels ++;
						}
						iNrCheckedPixels ++;
					}
				}
				iNewPoints += iIterationNewPoints;
				iIteration ++;

				//Dispose of borderregion as we clone it on every iteration
				oBorderRegion.Dispose();

				Debug.WriteLine(" Checked " + iNrCheckedPixels + " pixels, "  + iNrBorderPixels + " on border, included " + iIterationNewPoints + " new pixels",  "GrowRegion");


				//Do window events!
				Application.DoEvents();
			} while ((iIterationNewPoints > 0) && ((iIteration < iMaxNrIterations) || (iMaxNrIterations == -1)));

			Debug.WriteLine("Added a total of " + iNewPoints + " to the region", "GrowRegion");

			//Cleanup
			g.Dispose();

			if (iIterationNewPoints > 0) 
			{
				return(false);
			}
			else
			{
				return(true);
			}
		}
	#endregion

	#region Mask image functions
		public static unsafe Bitmap ToMask(Region oRgn, int iWidth, int iHeight)
		{
			//Create a mask image from a region
			Bitmap objBitmap = new Bitmap(iWidth, iHeight, System.Drawing.Imaging.PixelFormat.Format1bppIndexed);
			System.Drawing.Imaging.BitmapData bitmapData = objBitmap.LockBits(new Rectangle(0,0,iWidth,iHeight), 
				System.Drawing.Imaging.ImageLockMode.ReadWrite,
				System.Drawing.Imaging.PixelFormat.Format1bppIndexed );
																					
			int iPtrPixel0 = bitmapData.Scan0.ToInt32();
			byte * byPixel;
			int x, y, iRemainderX;;
			for (y = 0; y < bitmapData.Height; y ++)
			{
				byPixel = (byte *)(iPtrPixel0 + y * bitmapData.Stride);
				iRemainderX = 0;
				for (x = 0; x < bitmapData.Width; x ++)
				{
					//Pixel is on in region
					if (oRgn.IsVisible(x,y))
					{
						SetBitPixel(byPixel,iRemainderX, true);
					}

					if (iRemainderX == 7)
					{
						iRemainderX = 0;
						byPixel ++;
					}
					else
					{
						iRemainderX ++;
					}
				}
			}
			objBitmap.UnlockBits(bitmapData);

			return(objBitmap);
		}

		public static unsafe Region GetRegion(Bitmap objBitmap)
		{
			//Create a region from a mask
			Region oRgn = new Region();
			oRgn.MakeEmpty();
			Rectangle oPtRect = new Rectangle(0,0,1,1);
			System.Drawing.Imaging.BitmapData bitmapData = objBitmap.LockBits(new Rectangle(0,0,objBitmap.Width,objBitmap.Height), 
				System.Drawing.Imaging.ImageLockMode.ReadOnly,
				System.Drawing.Imaging.PixelFormat.Format1bppIndexed );
																					
			int iPtrPixel0 = bitmapData.Scan0.ToInt32();
			byte * byPixel;
			int x, y, iRemainderX;;
			for (y = 0; y < bitmapData.Height; y ++)
			{
				byPixel = (byte *)(iPtrPixel0 + y * bitmapData.Stride);
				iRemainderX = 0;
				for (x = 0; x < bitmapData.Width; x ++)
				{
					//Pixel is on in region
					if (GetBitPixel(byPixel, iRemainderX))
					{
						oPtRect.X = x;
						oPtRect.Y = y;
						oRgn.Union(oPtRect);
					}

					if (iRemainderX == 7)
					{
						iRemainderX = 0;
						byPixel ++;
					}
					else
					{
						iRemainderX ++;
					}
				}
			}
			objBitmap.UnlockBits(bitmapData);

			return(oRgn);
		}

		public static unsafe ArrayList GetBorderPoints(Bitmap oMask)
		{
			//Create an arraylist of arrays of points from a mask. Orignal mask is not modified!
			//Based upon a paper by Fu Chang, University of Taipei in Taiwan
			ArrayList oBorderPts = new ArrayList();

			Rectangle oRect = new Rectangle(0,0,oMask.Width,oMask.Height);

			Bitmap oLabeledPoints = new Bitmap(oRect.Width, oRect.Height, PixelFormat.Format1bppIndexed);

			//This bitmap needs to be larger to track contours lying on the border of the image
			Bitmap oNoGoPoints = new Bitmap(oRect.Width + 2, oRect.Height + 2, PixelFormat.Format1bppIndexed);

			System.Drawing.Imaging.BitmapData oMaskData = oMask.LockBits(oRect, 
				ImageLockMode.ReadOnly, PixelFormat.Format1bppIndexed );
			System.Drawing.Imaging.BitmapData oLPData = oLabeledPoints.LockBits(oRect, 
				ImageLockMode.ReadOnly, PixelFormat.Format1bppIndexed );
			System.Drawing.Imaging.BitmapData oNGData = oNoGoPoints.LockBits(oRect, 
				ImageLockMode.ReadOnly, PixelFormat.Format1bppIndexed );
																					
			byte *byMaskPixel0 = (byte *)oMaskData.Scan0.ToPointer(), byMaskPixel;
			byte *byLPPixel0 = (byte *)oLPData.Scan0.ToPointer(), byLPPixel;
			byte *byNGPixel0 = (byte *)oNGData.Scan0.ToPointer();

			int x, y, iRemainderX;
			for (y = 0; y < oMaskData.Height; y ++)
			{
				byMaskPixel = byMaskPixel0 + y * oMaskData.Stride;
				byLPPixel = byLPPixel0 + y * oLPData.Stride;
				iRemainderX = 0;
				for (x = 0; x < oMaskData.Width; x ++)
				{
					//Pixel is on
					if (GetBitPixel(byMaskPixel, iRemainderX)) 
					{

						//Check if pixel is unlabeled and pixel above it is not on. If so we have a new external contour
						if ( (!GetBitPixel(byLPPixel, iRemainderX)) && 
							( (y == 0) || !(GetBitPixel(byMaskPixel - oMaskData.Stride, iRemainderX)) ))
						{
							//Track the external contour
							TrackContour(x, y, true, oMaskData, oLPData, oNGData, oBorderPts);
						}

						//Check if pixel below is unmarked and on (not the same as unlabeled!)
						//The current pixel maybe labeled or unlabeled
						//We need to offset x and y because the NoGo bitmap is 1 pixel larger
						if ((!(GetBitPixel(byNGPixel0, oNGData.Stride, x + 1, y + 2)) && 
								((y == (oMaskData.Height - 1)) || !(GetBitPixel(byMaskPixel + oMaskData.Stride, iRemainderX)))))
						{
							//We have an internal contour, track it.
							TrackContour(x, y, false, oMaskData, oLPData, oNGData, oBorderPts);
						}
					}

					if (iRemainderX == 7)
					{
						iRemainderX = 0;
						byMaskPixel ++;
						byLPPixel ++;
					}
					else
					{
						iRemainderX ++;
					}
				}
			}

			//Cleanup
			oMask.UnlockBits(oMaskData);
			oLabeledPoints.UnlockBits(oLPData);
			oNoGoPoints.UnlockBits(oNGData);

			return(oBorderPts);
		}

		private static unsafe void TrackContour(int x, int y, 
																						bool bExternalContour, 
																						BitmapData oMaskData, 
																						BitmapData oLPData,
																						BitmapData oNGData, 
																						ArrayList oBorderPts)
		{
			//Start tracking. The order of the search directions is VERY cleverly crafted and important.
			if (bExternalContour)
			{
				Debug.WriteLine("Starting external contour at " + x + "," + y, "TrackContour");
			}
			else
			{
				Debug.WriteLine("Starting internal contour at " + x + "," + y, "TrackContour");
			}
			ArrayList oSingleContourBorderPts = new ArrayList();
			Point oFirstPt = new Point(x,y), oCurrentPt = new Point(x,y), oNewPt = new Point(x,y);
			oSingleContourBorderPts.Add(new Point(x,y));
			int iSearchDirection, iStartSearchDirection;
			Point [] oSearchPts = {new Point(1,0), new Point(1,1), new Point(0,1), new Point(-1,1),
															new Point(-1,0), new Point(-1,-1), new Point(0,-1), new Point(1,-1)};

			byte *byMaskPixel0 = (byte *)oMaskData.Scan0.ToPointer();
			byte *byLPPixel0 = (byte *)oLPData.Scan0.ToPointer();
			byte *byNGPixel0 = (byte *)oNGData.Scan0.ToPointer();
			
			if (bExternalContour)
			{
				iStartSearchDirection = 7;
			}
			else
			{
				iStartSearchDirection = 3;
			}

			do
			{
				//Do clockwise search until we find an 'On' pixel
				iSearchDirection = iStartSearchDirection;
				bool bIsolatedPt = true;
				do
				{
					oNewPt.X = oCurrentPt.X + oSearchPts[iSearchDirection].X;
					oNewPt.Y = oCurrentPt.Y + oSearchPts[iSearchDirection].Y;

					//Make sure pixel is inside image
					if ((oNewPt.X >= 0) && (oNewPt.X < oMaskData.Width) && (oNewPt.Y >= 0) && (oNewPt.Y < oMaskData.Height ))
					{
						//Break if we find an on pixel
						if (GetBitPixel(byMaskPixel0, oMaskData.Stride, oNewPt.X, oNewPt.Y)) 
						{
							bIsolatedPt = false;
							break;
						}
						else
						{
							//Mark pixel, so we avoid tracking the same contour twice
							SetBitPixel(byNGPixel0, oNGData.Stride, oNewPt.X + 1, oNewPt.Y + 1, true);
						}
					}
					else
					{
						//Mark border pixels (NoGO is 2 pixels larger than the mask), so we avoid tracking the same contour twice
						SetBitPixel(byNGPixel0, oNGData.Stride, oNewPt.X + 1, oNewPt.Y + 1, true);
					}

					iSearchDirection = (iSearchDirection + 1) % 8;

				}while (iSearchDirection != iStartSearchDirection);

				//We either have an isolated point or a new contour point
				if (bIsolatedPt)
				{
					//If the contour already contains points this is an error!
					if (oSingleContourBorderPts.Count > 1)
					{
						Debug.WriteLine("Found isolated point in non-singleton contour!", "TrackContour");
					}
					//Jump out of contour finding loop
					Debug.WriteLine("Stopped search at " + x + "," + y, "TrackContour");
					break;
				}
				else
				{
					//Add the point to the chain, and set all necessary variables for searching the next point
					//Also set pixel in the bitmaks of already labeled pixels, and its off neighbour in 
					//the NoGo image (this avoids tracking the same contour twice)
					//Debug.WriteLine("Added point " + oNewPt.X + "," + oNewPt.Y, "TrackContour");
					SetBitPixel(byLPPixel0, oLPData.Stride, oNewPt.X, oNewPt.Y, true);
					oCurrentPt.X = oNewPt.X;
					oCurrentPt.Y = oNewPt.Y;
					oSingleContourBorderPts.Add(new Point(oNewPt.X,oNewPt.Y));
					iStartSearchDirection = (iSearchDirection + 6) % 8;
				}
			}
			while ((oCurrentPt.X != oFirstPt.X) || (oCurrentPt.Y != oFirstPt.Y));

			//We now have a closed contour, we add it to the array of points collection
			Point [] oPts = new Point[oSingleContourBorderPts.Count];
			oSingleContourBorderPts.CopyTo(oPts);
			oBorderPts.Add(oPts);

			Debug.WriteLine("Added " + oSingleContourBorderPts.Count + "  points", "TrackContour");
		}

		public static unsafe void Border(Bitmap objBitmap)
		{
			Rectangle oBitmapRect = new Rectangle(0,0,objBitmap.Width, objBitmap.Height);
			System.Drawing.Imaging.BitmapData bitmapData = objBitmap.LockBits(oBitmapRect, 
				System.Drawing.Imaging.ImageLockMode.ReadWrite,
				System.Drawing.Imaging.PixelFormat.Format1bppIndexed);

			Border(bitmapData, oBitmapRect);
			objBitmap.UnlockBits(bitmapData);
		}

		private static unsafe void Border(BitmapData bitmapData, Rectangle oRect)
		{
			//Remove any but border pixels from a Mask
			Debug.WriteLine("Find border pixels " + oRect.ToString(),"ImageProcessing.Border");

			byte *byPixel0 = (byte*) bitmapData.Scan0.ToPointer();

			//Allocate and copy memory
			int iNrBytes = bitmapData.Stride*bitmapData.Height;
			IntPtr iPtr = Marshal.AllocHGlobal(iNrBytes);
			byte *byOriginalPixel0 = (byte *) iPtr.ToPointer();
			CopyMemory(byOriginalPixel0, byPixel0, iNrBytes);

			byte *byPixel, byOriginalPixel;
			int x, y, iRemainderX, xMin, yMin, xMax, yMax;
			yMin = (int) Math.Max(1, oRect.Top);
			yMax = (int) Math.Min(oRect.Bottom, bitmapData.Height - 1);
			xMin = (int) Math.Max(oRect.Left, 1);
			xMax = (int) Math.Min(oRect.Right, bitmapData.Width - 1);

			for (y = yMin; y < yMax; y ++)
			{
				int iOffset = y * bitmapData.Stride + xMin / 8;
				byPixel = byPixel0 + iOffset;
				byOriginalPixel = byOriginalPixel0 + iOffset;
				iRemainderX = xMin % 8;
				for (x = xMin; x < xMax; x ++)
				{
					if (GetBitPixel(byOriginalPixel, iRemainderX))
					{
						//Check connectedness.Use random access mode for next and previous pixel
						//Too avoid a lot of exception handling when crossing byte boundaries
						if (GetBitPixel(byOriginalPixel - bitmapData.Stride, iRemainderX) &&
							GetBitPixel(byOriginalPixel + bitmapData.Stride, iRemainderX) &&
							GetBitPixel(byOriginalPixel0, bitmapData.Stride, x + 1, y) && 
							GetBitPixel(byOriginalPixel0, bitmapData.Stride, x - 1, y) )
						{
							//Pixel is not on edge (not fully connected)
							SetBitPixel(byPixel, iRemainderX, false);
						}
					}

					if (iRemainderX == 7)
					{
						iRemainderX = 0;
						byPixel ++;
						byOriginalPixel ++;
					}
					else
					{
						iRemainderX ++;
					}
				}
			}
			Marshal.FreeHGlobal(iPtr);
		}

		public static unsafe void Morphology(Bitmap objBitmap, BitmapMaskMorphologyOperation iOperation)
		{
			Morphology(objBitmap, iOperation, new  Rectangle(0,0,objBitmap.Width,objBitmap.Height), 1);
		}

		public static unsafe void Morphology(Bitmap objBitmap, BitmapMaskMorphologyOperation iOperation, int iNrRepeats)
		{
			Morphology(objBitmap, iOperation, new  Rectangle(0,0,objBitmap.Width,objBitmap.Height), iNrRepeats);
		}
		
		private static unsafe void Morphology(Bitmap objBitmap, BitmapMaskMorphologyOperation iOperation, Rectangle oRect, int iNrRepeats)
			{
			//Erode/dilate a certain bitplane. This works MUCH faster than with a region, and we can still 
			//create a region from it afterwards. 
			//We always lock the whole bitmapData, and NOT the area defined by the rectangle ...
			int i;
			Rectangle oBitmapRect = new Rectangle(0,0,objBitmap.Width, objBitmap.Height);
			System.Drawing.Imaging.BitmapData bitmapData = objBitmap.LockBits(oBitmapRect, 
				System.Drawing.Imaging.ImageLockMode.ReadWrite,
				System.Drawing.Imaging.PixelFormat.Format1bppIndexed);


			switch(iOperation)
			{
				case BitmapMaskMorphologyOperation.Dilate:
					for (i=0; i < iNrRepeats; i ++) Morphology(bitmapData, false, oRect);
					break;
				case BitmapMaskMorphologyOperation.Erode:
					for (i=0; i < iNrRepeats; i ++) Morphology(bitmapData, true, oRect);
					break;
				case BitmapMaskMorphologyOperation.Open :
					for (i=0; i < iNrRepeats; i ++) Morphology(bitmapData, true, oRect);
					for (i=0; i < iNrRepeats; i ++) Morphology(bitmapData, false, oRect);
					break;
				case BitmapMaskMorphologyOperation.Close:
					for (i=0; i < iNrRepeats; i ++) Morphology(bitmapData, false, oRect);
					for (i=0; i < iNrRepeats; i ++) Morphology(bitmapData, true, oRect);
					break;
			}
			objBitmap.UnlockBits(bitmapData);
		}

		public static unsafe void Copy(Bitmap srcBitmap, Bitmap dstBitmap)
		{
			//Quick copying of data from one bitmap to another existing bitmap. They MUST be of 
			//the same type and size!!!!!!
			Rectangle oBitmapRect = new Rectangle(0,0,srcBitmap.Width, srcBitmap.Height);
			System.Drawing.Imaging.BitmapData srcBitmapData = srcBitmap.LockBits(oBitmapRect, 
				System.Drawing.Imaging.ImageLockMode.ReadOnly, srcBitmap.PixelFormat);
			System.Drawing.Imaging.BitmapData dstBitmapData = dstBitmap.LockBits(oBitmapRect, 
				System.Drawing.Imaging.ImageLockMode.ReadOnly, srcBitmap.PixelFormat);
			int nrBytes = srcBitmapData.Height * srcBitmapData.Stride;

			if (nrBytes == (dstBitmapData.Height * dstBitmapData.Stride) )
			{
				CopyMemory((byte *)dstBitmapData.Scan0.ToPointer(), (byte *) srcBitmapData.Scan0.ToPointer(), nrBytes);
			}
			else
			{
				MessageBox.Show(null, "Error copying bitmap data", "ImageProcessing.Copy", MessageBoxButtons.OK, MessageBoxIcon.Error);
			}

			//Release
			srcBitmap.UnlockBits(srcBitmapData);
			dstBitmap.UnlockBits(dstBitmapData);
		}

		private static unsafe void Morphology(BitmapData bitmapData, bool bErode, Rectangle oRect)
		{
			//Erode/dilate a certain bitplane. This works MUCH faster than with a region, and we can still 
			//create a region from it afterwards. This version can be used inside other functions which use
			//the bitmapdata, thereby avoiding lockouts and repeated marshaling of bitmap data
			//The rectangle can be used to further restrict the rectangle inside the bitmapdata over which 
			//iteration takes place, e.g. when being called from a function which already uses bitmapdata
			//over a larger rectangle.
			//We always assume the whole bitmap has been locked, and NOT only data in the rectangle
			Debug.WriteLine("Mask morphology in " + oRect.ToString(),"ImageProcessing.Morphology");

			byte *byPixel0 = (byte*) bitmapData.Scan0.ToPointer();

			//Allocate and copy memory
			int iNrBytes = bitmapData.Stride*bitmapData.Height;
			IntPtr iPtr = Marshal.AllocHGlobal(iNrBytes);
			byte *byOriginalPixel0 = (byte *) iPtr.ToPointer();
			CopyMemory(byOriginalPixel0, byPixel0, iNrBytes);

			byte *byPixel, byOriginalPixel;
			int x, y, iRemainderX, xMin, yMin, xMax, yMax;
			yMin = (int) Math.Max(1, oRect.Top);
			yMax = (int) Math.Min(oRect.Bottom, bitmapData.Height - 1);
			xMin = (int) Math.Max(oRect.Left, 1);
			xMax = (int) Math.Min(oRect.Right, bitmapData.Width - 1);

			for (y = yMin; y < yMax; y ++)
			{
				int iOffset = y * bitmapData.Stride + xMin / 8;
				byPixel = byPixel0 + iOffset;
				byOriginalPixel = byOriginalPixel0 + iOffset;
				iRemainderX = xMin % 8;
				for (x = xMin; x < xMax; x ++)
				{
					//Pixel needs to be on for erosion, off for dilation
					if (GetBitPixel(byOriginalPixel, iRemainderX) == bErode)
					{
						//Check connectedness.Use random access mode for next and previous pixel
						//Too avoid a lot of exception handling when crossing byte boundaries
						if (GetBitPixel(byOriginalPixel - bitmapData.Stride, iRemainderX)!=bErode || 
							GetBitPixel(byOriginalPixel + bitmapData.Stride, iRemainderX)!=bErode || 
							GetBitPixel(byOriginalPixel0, bitmapData.Stride, x + 1, y)!=bErode || 
							GetBitPixel(byOriginalPixel0, bitmapData.Stride, x - 1, y)!=bErode )
						{
							//Pixel is on edge (not fully connected), remove it for erosion, add it for dilation
							SetBitPixel(byPixel, iRemainderX, !bErode);
						}
					}

					if (iRemainderX == 7)
					{
						iRemainderX = 0;
						byPixel ++;
						byOriginalPixel ++;
					}
					else
					{
						iRemainderX ++;
					}
				}
			}
			Marshal.FreeHGlobal(iPtr);
		}

		public static unsafe void Invert(Bitmap mask)
		{
			//Convenience function, this is actually a Complement combine operation with an empty mask
			Bitmap tempMask = new Bitmap(mask.Width, mask.Height, PixelFormat.Format1bppIndexed);
			Combine(mask, tempMask, BitmapMaskCombineOperation.Complement);
			//mask.Save("C:\\temp\\InvertedMask.png");
			tempMask.Dispose();
		}

		public static unsafe void Combine(Bitmap oBitmap1, Bitmap oBitmap2, BitmapMaskCombineOperation iOperation)
		{
			Rectangle oRect = new Rectangle(0,0,oBitmap1.Width,oBitmap1.Height);
			Combine(oBitmap1, oBitmap2, iOperation, oRect);
		}

		public static unsafe void Combine(Bitmap oBitmap1, Bitmap oBitmap2, BitmapMaskCombineOperation iOperation, Rectangle oRect)
		{
			Rectangle oBitmapBB = new  Rectangle(0,0,oBitmap1.Width,oBitmap1.Height);
			System.Drawing.Imaging.BitmapData bitmap2Data = oBitmap2.LockBits(oBitmapBB, 
				System.Drawing.Imaging.ImageLockMode.ReadOnly ,
				System.Drawing.Imaging.PixelFormat.Format1bppIndexed);
			System.Drawing.Imaging.BitmapData bitmap1Data = oBitmap1.LockBits(oBitmapBB, 
				System.Drawing.Imaging.ImageLockMode.ReadWrite,
				System.Drawing.Imaging.PixelFormat.Format1bppIndexed);

			Combine(bitmap1Data, bitmap2Data, iOperation, oRect);
																		
			oBitmap1.UnlockBits(bitmap1Data);
			oBitmap2.UnlockBits(bitmap2Data);
		}

		public static unsafe void Combine(BitmapData bitmap1Data, BitmapData bitmap2Data, BitmapMaskCombineOperation iOperation,Rectangle oRect)
		{
			//This version can be used from other routines were the bitmap data has already been obtained																					
			byte * byPixel1_0 = (byte *)bitmap1Data.Scan0.ToPointer(),byPixel2_0 =  (byte *)bitmap2Data.Scan0.ToPointer();
			byte * byPixel1, byPixel2;
			int x, y, iRemainderX, xMin, yMin, xMax, yMax;

			yMin = (int) Math.Max(0, oRect.Top);
			yMax = (int) Math.Min(oRect.Bottom, bitmap1Data.Height);
			xMin = (int) Math.Max(oRect.Left, 0);
			xMax = (int) Math.Min(oRect.Right, bitmap1Data.Width);

			Debug.WriteLine("Combining masks in " + oRect.ToString(),"ImageProcessing.Combine");
			for (y = yMin; y < yMax; y ++)
			{
				int iOffset = y * bitmap1Data.Stride + xMin / 8;

				byPixel1 = byPixel1_0 + iOffset;
				byPixel2 = byPixel2_0 + iOffset;
				iRemainderX = xMin % 8;
				for (x = xMin; x < xMax; x ++)
				{
					bool bSet = false;
					switch(iOperation)
					{
						case BitmapMaskCombineOperation.Union:
							bSet = (GetBitPixel(byPixel1, iRemainderX) || GetBitPixel(byPixel2, iRemainderX));
							break;
						case  BitmapMaskCombineOperation.Intersect: 
							bSet = (GetBitPixel(byPixel1, iRemainderX) && GetBitPixel(byPixel2, iRemainderX));
							break;
						case   BitmapMaskCombineOperation.Xor:
							bool p1 = GetBitPixel(byPixel1, iRemainderX),p2 = GetBitPixel(byPixel2, iRemainderX);
							bSet = ((p1 && !p2) || (!p1 && p2));
							break;
						case   BitmapMaskCombineOperation.Exclude:
							bSet = (GetBitPixel(byPixel1, iRemainderX) && !GetBitPixel(byPixel2, iRemainderX));
							break;
						case   BitmapMaskCombineOperation.Complement:
							bSet = (!GetBitPixel(byPixel1, iRemainderX) && !GetBitPixel(byPixel2, iRemainderX));
							break;
					}
					SetBitPixel(byPixel1, iRemainderX, bSet);

					if (iRemainderX == 7)
					{
						iRemainderX = 0;
						byPixel1 ++;
						byPixel2 ++;
					}
					else
					{
						iRemainderX ++;
					}
				}
			}
		}

		public static unsafe bool GetBitPixel(byte * byBitplane, int iStride, int x, int y)
		{
			return (GetBitPixel(byBitplane + y * iStride + x / 8, x % 8));
		}

		public static unsafe bool GetBitPixel(byte * byBitPlane, int iRemainderX)
		{
			if ((*byBitPlane & (128 >> iRemainderX)) > 0)
			{
				return (true);
			}
			else
			{
				return (false);
			}
		}

		public static unsafe void SetBitPixel(byte * byBitPlane, int iStride, int x, int y, bool bOn)
		{
			SetBitPixel(byBitPlane + y * iStride + x / 8, x % 8, bOn);
		}

		public static unsafe void SetBitPixel(byte * byBitPlane, int iRemainderX, bool bOn)
		{
			byte iBitPlaneMask = (byte) (128 >> iRemainderX); //Other ordering, the one GDI + uses
			if (bOn)
			{
				*byBitPlane = (byte) (*byBitPlane | iBitPlaneMask);
			}
			else
			{
				*byBitPlane = (byte) (*byBitPlane & ~iBitPlaneMask);
			}
		}

		public unsafe static bool Grow(Bitmap oBitmap, 
																	 Bitmap oMask, 
																	 Point oSeedPt, 
																	 float sdEThreshold,
																	 bool useChromaticInfoOnly,
																	 int iMaxNrIterations)
		{
			//Grow a Mask using given ColorUtil and dE threshold given.
			//There is not much point in using unmanaged code as pixel access will be random
			//Returns after given number of iterations, with true when finished, false otherwise
			//This allows update of UI.
			//REMARK: detect boundingbox on first run, and use to limit 
			//number of iterations on subsequent runs ...
			Rectangle oBitmapBB = new  Rectangle(0,0,oBitmap.Width,oBitmap.Height);
			BitmapData maskData = oMask.LockBits(oBitmapBB, 
				System.Drawing.Imaging.ImageLockMode.ReadWrite,
				System.Drawing.Imaging.PixelFormat.Format1bppIndexed);
			byte * byMaskPtrPixel0 = (byte *) maskData.Scan0.ToPointer();
			Color oColor;
			ColorUtil oColorUtil = new ColorUtil();
			float[] sRGB = new float[3];
			float[] sLab = new float[3];
			float [] sSeedLab = new float[3];
			float fdE;

			//Compute ColorUtil of seed point in 3 by 3 area
			Rectangle oAvgRect = new Rectangle(oSeedPt.X - 1, oSeedPt.Y - 1, 3, 3);
			Average(oBitmap, oAvgRect, sRGB);
			ColorUtil.GammaRGBToRGB(sRGB, sRGB);
	    oColorUtil.sRGBToCIELab(sRGB, sSeedLab);
	    Debug.WriteLine("Seed point at " + oSeedPt.X + "," + oSeedPt.Y + " with CIE L*a*b* " + ColorUtil.ToString(sSeedLab));

			//Make sure seed point is in mask
			SetBitPixel(byMaskPtrPixel0, maskData.Stride, oSeedPt.X, oSeedPt.Y, true);

			int iIteration = 0, iNewPoints = 0, iIterationNewPoints;
			Rectangle oCurrentBB = new Rectangle(0,0,oBitmapBB.Width, oBitmapBB.Height);
	
			//The border pixels to be inspected.
			Bitmap oBorderMask;
			oBorderMask = (Bitmap) oMask.Clone();
			BitmapData maskBorderData = oBorderMask.LockBits(oBitmapBB, 
				System.Drawing.Imaging.ImageLockMode.ReadWrite,
				System.Drawing.Imaging.PixelFormat.Format1bppIndexed);
			byte * byBorderMaskPtrPixel0 = (byte *)maskBorderData.Scan0.ToPointer();

			//The pixels of the current region, plus those that have already been rejected!
			Bitmap oExclude;
			oExclude = (Bitmap) oMask.Clone();
			BitmapData ExcludeData = oExclude.LockBits(oBitmapBB, 
				System.Drawing.Imaging.ImageLockMode.ReadWrite,
				System.Drawing.Imaging.PixelFormat.Format1bppIndexed);
			byte * byExcludePtrPixel0 = (byte *)ExcludeData.Scan0.ToPointer();

			do
			{
				//Determine exterior border of region.
				Morphology(maskBorderData, false, oCurrentBB);
				Combine(maskBorderData, ExcludeData, BitmapMaskCombineOperation.Exclude, oCurrentBB);
				
				//Go through all the pixels, and add them to the original region if their ColorUtil 
				//is close enough to that of the seed point
 
				int x, y, iNrBorderPixels = 0, iNrCheckedPixels = 0, iRemainderX, yNewMin, yNewMax, xNewMax, xNewMin;
				byte * byBorderMaskPixel;
 
				yNewMin = (int) oCurrentBB.Height;
				yNewMax = (int) 0;
				xNewMax = (int) 0;
				xNewMin = (int) oCurrentBB.Width;
				iIterationNewPoints = 0;
				for (y = oCurrentBB.Top ; y < oCurrentBB.Bottom; y ++)
				{
					byBorderMaskPixel = byBorderMaskPtrPixel0 + y * maskBorderData.Stride + oCurrentBB.Left / 8;
					iRemainderX = oCurrentBB.Left % 8;
					for (x = oCurrentBB.Left; x < oCurrentBB.Right; x ++)
					{
						if (GetBitPixel(byBorderMaskPixel, iRemainderX))
						{
							//Keep track of bounding box. We have used separate variables because we cannot intialize a Rect to any values we choose
							if (x > xNewMax) xNewMax = x;
							if (x < xNewMin) xNewMin = x;
							if (y > yNewMax) yNewMax = y;
							if (y < yNewMin) yNewMin = y;

							//Check ColorUtil
							oColor = oBitmap.GetPixel(x,y);
							ColorUtil.ColorToRGB(oColor, sRGB);
							oColorUtil.sRGBToCIELab(sRGB, sLab);
		
							if (useChromaticInfoOnly)
							{
								fdE = ColorUtil.CIELabTodC(sLab, sSeedLab);
							}
							else
							{
								fdE = ColorUtil.CIELabTodE(sLab, sSeedLab);
							}
							//Debug.WriteLine("Pixel at " + x + "," + y + " has dE " + fdE + " with seed", "GrowRegion");

							//Add to region if below threshold
							if (fdE < sdEThreshold)
							{
								SetBitPixel(byMaskPtrPixel0, maskData.Stride, x, y, true);
								iIterationNewPoints ++;
							}
							else
							{
								SetBitPixel(byBorderMaskPixel, iRemainderX, false);
							}

							//Add to exclusion region
							SetBitPixel(byExcludePtrPixel0, ExcludeData.Stride, x, y, true);

							//Count number of border pixels
							iNrBorderPixels ++;
						}

						if (iRemainderX == 7)
						{
							iRemainderX = 0;
							byBorderMaskPixel ++;
						}
						else
						{
							iRemainderX ++;
						}

						//Count number of pixels checked
						iNrCheckedPixels ++;
					}
				}

				//Update the boundingbox
				oCurrentBB.Y = Math.Max(yNewMin - 2,0);
				oCurrentBB.Height  = Math.Min(yNewMax + 2, oBitmapBB.Bottom) - oCurrentBB.Y;
				oCurrentBB.X = Math.Max(xNewMin - 2,0);
				oCurrentBB.Width  = Math.Min(xNewMax + 2, oBitmapBB.Right) - oCurrentBB.X;

				iNewPoints += iIterationNewPoints;
				iIteration ++;

				Debug.WriteLine(" Checked " + iNrCheckedPixels + " pixels, "  + iNrBorderPixels + " on border, included " + iIterationNewPoints + " new pixels",  "GrowRegion");

				//Do window events
				//Application.DoEvents();
			} while ((iIterationNewPoints > 0) && ((iIteration < iMaxNrIterations) || (iMaxNrIterations == -1)));

			Debug.WriteLine("Added a total of " + iNewPoints + " to the region", "GrowRegion");
			oBorderMask.UnlockBits(maskBorderData);
			oMask.UnlockBits(maskData);
			oExclude.UnlockBits(ExcludeData);

			if (iIterationNewPoints > 0) 
			{
				return(false);
			}
			else
			{
				return(true);
			}
		}

		public unsafe static void Find(Bitmap oBitmap, 
																	 Bitmap oMask, 
																	 float[] sSeedLab, 
																	 float sdEThreshold)
		{
			//Find all pixels within a dE rangle from given Lab ColorUtil
			Rectangle oRect = new  Rectangle(0,0,oBitmap.Width,oBitmap.Height);
			BitmapData maskData = oMask.LockBits(oRect, 
				System.Drawing.Imaging.ImageLockMode.ReadWrite,
				System.Drawing.Imaging.PixelFormat.Format1bppIndexed);
			byte * byMaskPtrPixel0 = (byte *) maskData.Scan0.ToPointer();

			System.Drawing.Imaging.BitmapData bitmapData = oBitmap.LockBits(oRect, 
				System.Drawing.Imaging.ImageLockMode.ReadOnly ,
				System.Drawing.Imaging.PixelFormat.Format24bppRgb);
			byte * byBitmapPtrPixel0 = (byte *) bitmapData.Scan0.ToPointer(), byBitmapPtrPixel;
																					
			ColorUtil oColorUtil = new ColorUtil();
			int[] iRGB = new int[3];
			float[] sRGB = new float[3];
			float[] sLab = new float[3];
			float fdE;

				//Go through all the pixels, and add them to the original region if their ColorUtil 
				//is close enough to that of the given ColorUtil
				int x, y;
				for (y = 0; y < oBitmap.Height; y ++)
				{
					byBitmapPtrPixel = byBitmapPtrPixel0 + y * bitmapData.Stride;

					for (x = 0; x < oBitmap.Width; x ++)
					{

					//Check ColorUtil
					iRGB[2] = *byBitmapPtrPixel;
					byBitmapPtrPixel ++;
					iRGB[1] = *byBitmapPtrPixel;
					byBitmapPtrPixel ++;
					iRGB[0] = *byBitmapPtrPixel;
					byBitmapPtrPixel ++;

					ColorUtil.GammaRGBToRGB(iRGB, sRGB);
					oColorUtil.sRGBToCIELab(sRGB, sLab);
					fdE = ColorUtil.CIELabTodE(sLab, sSeedLab);
						
					//Add to region if below threshold
					if (fdE < sdEThreshold)
					{
						SetBitPixel(byMaskPtrPixel0, maskData.Stride, x, y, true);
					}
				}
			}

			oMask.UnlockBits(maskData);
			oBitmap.UnlockBits(bitmapData);
		}

		#endregion

	}
}
